For this project we will be looking at LiDAR-generated Canopy Height Models and looking at uncertainty compared to field measurements of plot-level and tree-level characteristics.
For this project, I am going to be working with data from NEON Domain 01, Harvard Forest Field Site (HARV). I downloaded my data from NEON’s figshare site.
Specific characteristics we will be looking at include… *Fill this part in
Also, added additional setup steps in this code chunk.
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# also set stringsasFactors to false
options(stringsAsFactors = FALSE)
site_chm <- raster("../NEONdata/D01-Massachusetts/HARV/2014/lidar/HARV_lidarCHM.tif")
# take a look at it
plot(site_chm, main="Canopy Height Model\nHARV Field Site")
hist(site_chm, main="CHM Histogram\nHARV Field Site")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 10% of the raster cells were used. 100000 values used.
# set 0 values to NA
site_chm[site_chm == 0] <- NA
# take another look at the histogram
hist(site_chm, main="CHM Histogram\nHARV Field Site",
xlab="Height (m)")
site_insitu <- read.csv("../NEONdata/D01-Massachusetts/HARV/2014/insitu/veg_structure/D01_2012_HARV_vegStr.csv",
header=TRUE, sep=",")
# let's make sure the data came in okay
head(site_insitu) # lots of NA values
## site_id sitename plotid easting northing taxonid
## 1 HARV Harvard Forest A01 NA NA NA
## 2 HARV Harvard Forest A01 NA NA NA
## 3 HARV Harvard Forest A01 NA NA NA
## 4 HARV Harvard Forest A01 NA NA NA
## 5 HARV Harvard Forest A01 NA NA NA
## 6 HARV Harvard Forest A01 NA NA NA
## scientificname individual_id pointid individualdistance
## 1 Vaccinium corymbosum NA 41 NA
## 2 Vaccinium corymbosum NA 42 NA
## 3 Acer rubrum NA 43 NA
## 4 Castanea dentata NA 44 NA
## 5 Castanea dentata NA 45 NA
## 6 Tsuga canadensis NA 46 NA
## individualazimuth dbh dbhheight basalcanopydiam basalcanopydiam_90deg
## 1 NA 0.5 NA NA NA
## 2 NA 0.4 NA NA NA
## 3 NA 5.3 NA NA NA
## 4 NA 6.5 NA NA NA
## 5 NA 8.3 NA NA NA
## 6 NA 7.5 NA NA NA
## maxcanopydiam canopydiam_90deg canopy_diameter vertical_position
## 1 NA 0.0 0.7 Understory
## 2 NA 0.0 0.8 Understory
## 3 NA 5.4 3.7 Understory
## 4 NA 0.0 5.0 Understory
## 5 NA 7.7 0.0 Understory
## 6 NA 0.0 0.0 Understory
## stem_height stemremarks stemstatus canopyform livingcanopy inplotcanopy
## 1 1.5 NA NA NA NA NA
## 2 1.6 NA NA NA NA NA
## 3 10.5 NA NA NA NA NA
## 4 7.0 NA NA NA NA NA
## 5 2.4 NA NA NA NA NA
## 6 4.4 NA NA NA NA NA
## materialsampleid dbhqf stemmapqf plant_group common_name aop_plot
## 1 NA NA NA NA NA NA
## 2 NA NA NA NA NA NA
## 3 NA NA NA NA NA NA
## 4 NA NA NA NA NA NA
## 5 NA NA NA NA NA NA
## 6 NA NA NA NA NA NA
## unique_id
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA
# look at a histogram of the heights
hist(site_insitu$stem_height, main="In Situ Histogram\nHARV Field Site",
xlab="Height (m)")
In the insitu data, it looks like there’s a lot more shorter vegetation than was captured in the CHM.
# loading using readOGR
site_plots <- readOGR("../NEONdata/D01-Massachusetts/HARV/vector_data/",
"HARV_PlotCentroids")
## OGR data source with driver: ESRI Shapefile
## Source: "../NEONdata/D01-Massachusetts/HARV/vector_data/", layer: "HARV_PlotCentroids"
## with 21 features
## It has 25 fields
# look at the plot locations
plot(site_plots)
# overlay plot locations on CHM
plot(site_chm, main="CHM with In Situ Locations\nHARV Field Site")
plot(site_plots, add=TRUE)
Add link here
I’m going to set up square plots around the plot centroid. This process involves setting a radius, which will be half the length of one side of the square.
This is based on the assumption that plots are oriented north and not rotated.
# first set radius
radius <- 20 # in meters
# then define the plot boundaries using the radius
yPlus <- site_plots$Y + radius
xPlus <- site_plots$X + radius
yMinus <- site_plots$Y - radius
xMinus <- site_plots$X - radius
ID <- site_plots$plotID
# let's look at our plot IDs
ID
## [1] "HARV_015" "HARV_033" "HARV_034" "HARV_035" "HARV_036" "HARV_037"
## [7] "HARV_038" "HARV_039" "HARV_040" "HARV_041" "HARV_042" "HARV_043"
## [13] "HARV_044" "HARV_045" "HARV_046" "HARV_047" "HARV_048" "HARV_049"
## [19] "HARV_050" "HARV_051" "HARV_052"
# here I am just making a list of the coordinates of each corner of the square
# will follow x, y, x, y, etc.
# later we will turn this into the 4 corners of the square the bounds the plot
# this requires 5 corners to draw the square - why are there 6?
square <- cbind(xMinus, yPlus,
xPlus, yPlus,
xPlus, yMinus,
xMinus, yMinus,
xMinus, yPlus,
xMinus, yPlus)
# Create a function to do this
polys <- SpatialPolygons(mapply(function(poly, id) {
xy <- matrix(poly, ncol=2, byrow=TRUE) # take a list and create a matrix
Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID),proj4string=CRS(as.character("+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")))
This is an important step to extract the data.
polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))
plot(site_chm, main="CHM with square in situ plots\nHARV Field Site")
plot(polys.df, add=TRUE)
# extracting max height from CHM at field data sites
chm_maxheight <- extract(site_chm, polys.df,
fun=max, sp=TRUE, # sp for return spatial object
stringsasFactors=FALSE)
# extract mean height from CHM at field data sites
chm_meanheight <- extract(site_chm, polys.df,
fun=mean, sp=TRUE, # sp for return spatial object
stringsasFactors=FALSE)
# check outputs
chm_maxheight@data
## id HARV_lidarCHM
## HARV_015 HARV_015 19.73
## HARV_033 HARV_033 30.06
## HARV_034 HARV_034 29.86
## HARV_035 HARV_035 23.20
## HARV_036 HARV_036 29.09
## HARV_037 HARV_037 25.58
## HARV_038 HARV_038 28.87
## HARV_039 HARV_039 25.74
## HARV_040 HARV_040 NA
## HARV_041 HARV_041 27.34
## HARV_042 HARV_042 26.47
## HARV_043 HARV_043 26.24
## HARV_044 HARV_044 23.92
## HARV_045 HARV_045 25.82
## HARV_046 HARV_046 20.02
## HARV_047 HARV_047 22.85
## HARV_048 HARV_048 26.85
## HARV_049 HARV_049 23.51
## HARV_050 HARV_050 28.59
## HARV_051 HARV_051 24.97
## HARV_052 HARV_052 28.87
chm_meanheight@data
## id HARV_lidarCHM
## HARV_015 HARV_015 15.76019
## HARV_033 HARV_033 24.50129
## HARV_034 HARV_034 22.06391
## HARV_035 HARV_035 16.70113
## HARV_036 HARV_036 24.13921
## HARV_037 HARV_037 19.59647
## HARV_038 HARV_038 21.68227
## HARV_039 HARV_039 21.75469
## HARV_040 HARV_040 NA
## HARV_041 HARV_041 19.79854
## HARV_042 HARV_042 21.94006
## HARV_043 HARV_043 19.49369
## HARV_044 HARV_044 17.43806
## HARV_045 HARV_045 21.03539
## HARV_046 HARV_046 12.43345
## HARV_047 HARV_047 15.44984
## HARV_048 HARV_048 19.45156
## HARV_049 HARV_049 18.45651
## HARV_050 HARV_050 20.35642
## HARV_051 HARV_051 20.80787
## HARV_052 HARV_052 19.81688
unique(site_insitu$plotid)
## [1] "A01" "A04" "A05" "A06" "A07" "A08" "A09" "A12" "A13" "A16" "A17"
## [12] "A18" "A19" "A20" "A31" "A34" "A10" "A15"
# extract max ht by plot from insitu data with dplyr
insitu_maxheight <- site_insitu %>%
group_by(plotid) %>%
summarise(max(stem_height))
names(insitu_maxheight) <- c("plotid","max_height")
hist(insitu_maxheight$max_height)
hist(chm_maxheight@data$HARV_lidarCHM)
# quick look at data summary
summary(fieldmaxht.sjer)
summary(height.sjer@data$SJER_lidarCHM)
# update names
names(fieldmaxht.sjer) <- c("plotid","field_maxht")